Librerias

library("rmdformats")
library("VGAM")
library("readr")
library("here")
library("Rsolnp")
library("GB2")
library("fitdistrplus")
library("tidyverse")
library("viridis")
library("gridExtra")

Practica con Pareto

Generamos datos random de la distribución, poniendo atención al tercer parámetro, que es el valor de alfa

dataset = rsinmad(10000, 2.07, 1.14, 1.75)

Sacamos los datos a graficar aplicando las operaciones pertinentes

log_de_x= log(dataset)
log_de_acum = log(1-psinmad(dataset, 2.07, 1.14, 1.75))

Hacemos los datasets para graficar la relación lineal utilizando el valor de alfa como la pendiente

val_x = seq(-4, 6, .1)
val_y = -1.75*(val_x -2) - 2.9961

Por últmo, graficamos esta relación lineal

plot(log_de_x, log_de_acum, main="ejemplo de diagrama de pareto")
lines(val_x, val_y, col="red")

Practica con Log normal

Primero se especifica la media y la mediana reportada por el U.S. Census Buerau. En este ejemplo se utiliza la del año 2008. La informaci?n puede encontrarse en https://www.census.gov/data/tables/time-series/demo/income-poverty/historical-income-households.html Las cantidades estan en miles de dolares.

media_oficial = 80
mediana_oficial = 50

Utilizando las relaciones entre la media y la varianza de la normal asociada calculamos los parametros para dicha normal

media = log(mediana_oficial)
media_mas_varianza = log(media_oficial)

Restamos la media a la media mas la varianza para poder obtener la varianza solita y procedemos a sacar la raiz para obtener la desviacion estandar

desv_est = sqrt(2*log(media_oficial/mediana_oficial))

Se hace un array con los cuantiles para los cuales se quiere estimar el ratio de desigualdad

quant = c(0, .1, .2, .3, .4, .5, .6, .7, .75, .8, .85, .9, .95, .98, .99, .995, .999, .9995, .9999, .99999)

Utilizando la media y la desviaci?n est?ndar que ya obtuvimos, podemos entonces calcular los ingresos por household que corresponden a cada cuantil

ing_por_household = qlnorm(quant, media, desv_est)

Posteriormente, utilizando la relaci?n entre la distribuci?n de hogares por ingreso y la distribuci?n del ingreso total, podemos calcular los par?metros de la distribuci?n del ingreso total de la siguiente forma

media_ing_total = media + (desv_est)^2
desv_est_total = desv_est

Con estos parámetros, podemos estimar el porcentaje del ingreso total por ingreso

ing_total_per = plnorm(ing_por_household, media_ing_total, desv_est_total)

Se calcula entonces el área que hay a partir de cada punto

quant_inverso = 1-quant
int_total_per_inverso = 1- ing_total_per

en quant_inverso se tiene el porcentaje de hogares que tienen un ingreso mayor a los puntos en ing_por_household y en int_total_per_inverso esta el porcentaje total del ingreso que acaparan los hogares que tienen un ingreso mayor a los puntos en ing_per_household

Para calcular el coeficiente relativo, unicamente dividimos el porcentaje del ingreso total acaparado por hogares con un ingreso mayor a un punto determinado entre el porcentaje de (el número) de casas que tienen un ingreso mayor al mismo punto.

coef_rel = int_total_per_inverso / quant_inverso

Por ultimo se pone todo en un dataframe para que tenga el formato de tabla

tabla_final = data.frame("cuantil" = quant,
                         "Ingreso en K dolares" = ing_por_household,
                         "Ingreso Total hasta k dolares" = ing_total_per,
                         "P(X>x)" = quant_inverso,
                         "P(Y>y)" = int_total_per_inverso,
                         "ratio"= coef_rel)

print(tabla_final)
##    cuantil Ingreso.en.K.dolares Ingreso.Total.hasta.k.dolares  P.X.x.
## 1  0.00000              0.00000                    0.00000000 1.00000
## 2  0.10000             14.43286                    0.01218988 0.90000
## 3  0.20000             22.11017                    0.03505797 0.80000
## 4  0.30000             30.07204                    0.06759562 0.70000
## 5  0.40000             39.11058                    0.11068622 0.60000
## 6  0.50000             50.00000                    0.16613799 0.50000
## 7  0.60000             63.92133                    0.23693621 0.40000
## 8  0.70000             83.13370                    0.32810958 0.30000
## 9  0.75000             96.15559                    0.38397786 0.25000
## 10 0.80000            113.07013                    0.44910674 0.20000
## 11 0.85000            136.57670                    0.52666683 0.15000
## 12 0.90000            173.21579                    0.62248424 0.10000
## 13 0.95000            246.35508                    0.75026183 0.05000
## 14 0.98000            366.21264                    0.86086397 0.02000
## 15 0.99000            476.99674                    0.91257891 0.01000
## 16 0.99500            607.52408                    0.94589485 0.00500
## 17 0.99900           1000.37074                    0.98302616 0.00100
## 18 0.99950           1214.78132                    0.98985623 0.00050
## 19 0.99990           1840.43491                    0.99701548 0.00010
## 20 0.99999           3124.42020                    0.99950851 0.00001
##          P.Y.y.     ratio
## 1  1.0000000000  1.000000
## 2  0.9878101217  1.097567
## 3  0.9649420294  1.206178
## 4  0.9324043786  1.332006
## 5  0.8893137759  1.482190
## 6  0.8338620125  1.667724
## 7  0.7630637877  1.907659
## 8  0.6718904199  2.239635
## 9  0.6160221442  2.464089
## 10 0.5508932573  2.754466
## 11 0.4733331743  3.155554
## 12 0.3775157585  3.775158
## 13 0.2497381734  4.994763
## 14 0.1391360276  6.956801
## 15 0.0874210871  8.742109
## 16 0.0541051484 10.821030
## 17 0.0169738402 16.973840
## 18 0.0101437723 20.287545
## 19 0.0029845247 29.845247
## 20 0.0004914939 49.149394

Simulación datos censurados

Simulamos primero datos aleatorios de una distribución Singh - Maddala, censurando artificialmente los datos que sean mayores a 3. De forma preeliminar, calculamos el porcentaje de los datos que estamos censurando y cual es la media real de los agentes con un ingreso mayor a 3

ingresos <- rsinmad(100000, 1.14, 1.75, 2.07)

top_code <- 3

prop_real <- sum(ingresos>top_code)/100000
real_mean <- mean(ingresos[ingresos>top_code])

Procedemos a hacer la censura de los datos

ingresos_censurado <- rep(0,100000)
ingresos_censurado[ingresos < top_code] <- ingresos[ingresos < top_code]
ingresos_censurado[ingresos >= top_code] <- top_code

Método

Suponemos entonces que la distribución de los ingresos para los datos arriba del cuantil 95% sigue una distribución Pareto Tipo I, y usando el estimador propuesto por Armour, Burkhauser y Larrimore estimamos el parámetro \(\alpha\).

Una vez que tenemos dicho parámetro, utilizamos la propiedad de la distribución pareto para calcular la media del ingreso de los agentes con un ingreso mayor a 3 y utilizamos esta media para imputar los datos que habían sido censurados. Esta nueva serie conservará mejor ciertas propiedades de los datos sin censura, como la media.

cutoff <- quantile(ingresos, .95)[[1]]
M = sum(ingresos_censurado>cutoff & ingresos_censurado<top_code)
Te = sum(ingresos_censurado==3)
alpha_est = M/(Te*log(top_code) + sum(log(ingresos_censurado[ingresos_censurado>=cutoff & ingresos_censurado<top_code])) - (M+Te)*log(cutoff))

invertido = alpha_est/(alpha_est-1)
estimated_mean <- top_code*invertido

## Media estimada
estimated_mean
## [1] 4.562036
## Media real
real_mean
## [1] 4.28969

Notamos que la media estimada y la real son muy similares.

Diferencia entre truncamiento y censura

x <- rnorm(100000,0,1)
y <- rnorm(100000,0,1)

t <- x[x<1.5]
c <- y
c[y>=1.5] <- 1.5

x <- x+1

ggplot()+
  geom_density(aes(x=t), fill = "#f04546", size = 1)+
  geom_density(aes(x=c), fill = "#6db823", size = 1) +
  geom_density(aes(x=x), fill = "#5b5c5a", size = 1, alpha=.4) +
  geom_vline(xintercept = 1.5, linetype = 2, color = "#3591d1", size = 1) +
  xlab("x") + ylab("Densidad")+
  theme_bw()

Practica con la ENIGH

Cargamos datos

## ENIGH
## Homogeneizamos los nombres de las columnas de interes
concentradohogar2012 <- read_csv(here("data", "concentradohogar2012.csv"))
colnames(concentradohogar2012)[colnames(concentradohogar2012) == "factor_hog"] <- "factor"

concentradohogar2014 <- read_csv(here("data", "concentradohogar2014.csv"))
colnames(concentradohogar2014)[colnames(concentradohogar2014) == "factor_hog"] <- "factor"

concentradohogar2016 <- read_csv(here("data", "concentradohogar2016.csv"))
concentradohogar2018 <- read_csv(here("data", "concentradohogar2018.csv"))

## Datos anonimos del SAT
sat2012 <- read_delim(here("data", "RUIDO_FINAL_PF_2012.tab"), "\t", escape_double = FALSE, trim_ws = TRUE)

Histogramas

Creamos histogramas con diferentes parámetros para familiarizarnos con los datos.

df <- tibble(ing_cor = c(concentradohogar2012$ing_cor,
                         concentradohogar2014$ing_cor,
                         concentradohogar2016$ing_cor,
                         concentradohogar2018$ing_cor),
             factor = c(concentradohogar2012$factor,
                        concentradohogar2014$factor,
                        concentradohogar2016$factor,
                        concentradohogar2018$factor),
o = c(rep("2012", length(concentradohogar2012$ing_cor)),
                     rep("2014", length(concentradohogar2014$ing_cor)),
                     rep("2016", length(concentradohogar2016$ing_cor)),
                     rep("2018", length(concentradohogar2018$ing_cor)))
              )

histograma <- ggplot(df, aes(x = ing_cor, weight = factor)) + 
  geom_histogram(bins = 60) + 
  facet_wrap(~año) +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::comma) +
  xlab("Ingreso Corriente") +  ylab("Frecuencia") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

histograma_zoomed <- ggplot(df, aes(x = ing_cor, weight = factor)) + 
  geom_histogram(bins = 600) + 
  facet_wrap(~año) +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::comma) +
  coord_cartesian(xlim=c(0,3000000)) +
  xlab("Ingreso Corriente") +  ylab("Frecuencia") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

histograma_log <- ggplot(df, aes(x = ing_cor, weight = factor)) + 
  geom_histogram(bins = 30) + 
  facet_wrap(~año) +
  scale_x_continuous(trans = "log", labels = scales::comma) +
  scale_y_continuous(labels = scales::comma) +
  xlab("Ingreso Corriente (log)") +  ylab("Frecuencia") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Histograma automático; parámetros usuales

histograma

Histograma con zoom; cortamos los ejes

histograma_zoomed

Histograma escala logaritmica

histograma_log

Ingreso medio por percentiles

Ya que los histogramas no nos fueron de gran utilidad, nos pasamos a investigar qué pasa en cada decil de la distribución de los ingresos. Además, buscamos ver si existe consistencia entre cada año de aplicación de la ENIGH.

ingreso_medio_percentil <- function(datos, n_percentil) {
  # datos son las base de datos del censo, n_percentil es el numero de divisiones (e.g.: 10 retorna deciles, 100 retorna percentiles)
  
  ## Datos sorteados por ingresos trimestales corrientes más indice de los hogares representados observados
  sorted_ing <- datos %>% select(ing_cor, factor) %>%
    arrange(ing_cor) %>%
    mutate(cumulative_obs = cumsum(factor))
  ## Total de hogares representados (ver definicion de factor en documentacion de la ENIGH)
  total_hogares <- sum(datos$factor)
  
  ## Tamaño por decil truncado para no incluir decimales
  tam_decile <- trunc(total_hogares/n_percentil)
  
  ## Primeros n-1 deciles
  lim_inf <- 1
  lim_sup <- tam_decile
  ing_med_de <- c()
  for (i in 1:(n_percentil-1)) {
    ing_med_de[i] <- sorted_ing %>% filter(cumulative_obs >= lim_inf, cumulative_obs <= lim_sup) %>%
      summarise(mean = sum(ing_cor*factor)/tam_decile) %>% 
      pull()
    
    ## Limites del decil
    lim_inf <- lim_inf + tam_decile
    lim_sup <- lim_sup + tam_decile
  }
  
  ## Decimo decil aparte para incluir el resto de observaciones
  ing_med_de[n_percentil] <- sorted_ing %>% filter(cumulative_obs >= lim_inf) %>%
    summarise(mean = sum(ing_cor*factor)/tam_decile) %>% 
    pull()
  
  return(ing_med_de)
}

n_percentil <- 10
df <- tibble("2012" = ingreso_medio_percentil(concentradohogar2012, n_percentil),
             "2014" = ingreso_medio_percentil(concentradohogar2014, n_percentil),
             "2016" = ingreso_medio_percentil(concentradohogar2016, n_percentil),
             "2018" = ingreso_medio_percentil(concentradohogar2018, n_percentil),
             index = 1:n_percentil)

df <- df %>% gather(key = "Año", value = "Ingreso_medio", -index)

ingreso_medio_decil <- ggplot(df , aes(x = index, y = Ingreso_medio))+
  geom_bar(stat="identity")+
  facet_wrap(~ Año) +
  scale_y_continuous(labels = scales::comma, breaks = seq(0,175000,25000))+
  scale_x_continuous(breaks = 1:10, labels = scales::ordinal)+
  xlab("Decil") + ylab("Ingreso medio") +
  labs(title="Ingreso medio por deciles") +
  theme_bw()

ingreso_medio_decil

Ingreso acumulado por Decil

Quizá nos sea más fácil observar el comportamiento por decil si lo graficamos en términos de los porcentajes que acumula cada decil. Además esto nos dará una idea muy clara de la desigualdad presente entre el decil más rico y el decil más pobre.

ingreso_acumulado_percentil <- function(datos, n_percentil){
  # Retorna el porcentaje del ingreso total que acumula cada percentil
  # datos son las base de datos del censo, n_percentil es el numero de divisiones (e.g.: 10 retorna deciles, 100 retorna percentiles)
  
  ## Datos sorteados por ingresos trimestales corrientes más indice de los hogares representados observados
  sorted_ing <- datos %>% select(ing_cor, factor) %>%
    arrange(ing_cor) %>%
    mutate(cumulative_obs = cumsum(factor))
  
  ## Total de hogares representados (ver definicion de factor en documentacion de la ENIGH)
  total_hogares <- sum(datos$factor)
  
  ## Muestras por por percentil truncado para no incluir decimales
  tam_decile <- trunc(total_hogares/n_percentil)
  
  
  ## Primeros n-1 deciles
  lim_inf <- 1
  lim_sup <- tam_decile
  ing_total <- c()
  for (i in 1:(n_percentil-1)) {
    ing_total[i] <- sorted_ing %>% filter(cumulative_obs >= lim_inf, cumulative_obs <= lim_sup) %>%
      summarise(mean = sum(ing_cor*factor)) %>% 
      pull()
    
    ## Limites del decil
    lim_inf <- lim_inf + tam_decile
    lim_sup <- lim_sup + tam_decile
  }
  
  ## Decimo decil aparte para incluir el resto de observaciones
  ing_total[n_percentil] <- sorted_ing %>% filter(cumulative_obs >= lim_inf) %>%
    summarise(mean = sum(ing_cor*factor)) %>% 
    pull()
  
  ing_total <- ing_total/sum(datos$ing_cor*datos$factor)
  
  return(ing_total)
}
n_percentil <- 10
df <- tibble("2012" = ingreso_acumulado_percentil(concentradohogar2012, n_percentil),
             "2014" = ingreso_acumulado_percentil(concentradohogar2014, n_percentil),
             "2016" = ingreso_acumulado_percentil(concentradohogar2016, n_percentil),
             "2018" = ingreso_acumulado_percentil(concentradohogar2018, n_percentil),
             index = fct_rev(as.factor(1:n_percentil))) 

df <- df %>% gather(key = "Año", value = "Ingreso", -index)

ingreso_acu_decil <- ggplot(df, aes(fill = index, x = Año, y = Ingreso)) +
  geom_bar(position="stack", stat="identity") +
  ylab("Porcentaje de Ingreso acumulado") + xlab("Año") +
  scale_fill_viridis(discrete = T, name = "Decil") +
  scale_y_continuous(labels = scales::percent) +
  theme_bw() +
  labs(title= "Ingreso acumulado por Decil")+
  theme(legend.background = element_rect(fill = "transparent", colour = "transparent"))

ingreso_acu_decil

Ingreso acumulado en los ultimos 0.1% de los hogares

Ya que esta ultima gráfica tampoco nos arroja información contundente sobre la consistencia de los datos a traves de los años, exploraremos qué ocurre con la acaparación de los ingresos en los ultimos cuantiles de los datos.

n_percentil <- 1000
df <- tibble("2012" = ingreso_acumulado_percentil(concentradohogar2012, n_percentil),
             "2014" = ingreso_acumulado_percentil(concentradohogar2014, n_percentil),
             "2016" = ingreso_acumulado_percentil(concentradohogar2016, n_percentil),
             "2018" = ingreso_acumulado_percentil(concentradohogar2018, n_percentil),
             index = 1:n_percentil)

df <- df %>% gather(key = "Año", value = "Ingreso", -index)

ingreso_acu_per <- ggplot(df, aes(x = index, y = Ingreso)) +
  geom_line(aes(color = Año), size = 1) +
  scale_x_continuous(labels = scales::ordinal) +
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(xlim=c(990,1000)) +
  ylab("Porcentaje de Ingreso acumulado") + xlab("Percentil (0.1%)") +
  labs(title="Ingreso acumulado por cada 0.1% de la poblacion") +
  theme_bw() +
  theme(legend.position = c(.1,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))

ingreso_acu_per

Gastan más de lo que ingresan

Otra medida que nos puede arrojar información sobre posibles incosistencias es saber si hay gente que reportó que gasta más de lo que dijo que ingresaba. De la misma manera, lo haremos para los 4 censos bajo estudio para saber si podemos confiar en la métrica dada.

gastan_de_mas <- function(datos){
  total_hogares <- sum(datos$factor)
  datos %>% 
    transmute(mas_gasto = (gasto_mon > ing_cor)*factor) %>%
    summarise(proporcion = sum(mas_gasto)/total_hogares) %>%
    pull()
}

# Porcentaje de personas que gastan de más
over_spenders <- c(gastan_de_mas(concentradohogar2012),
                   gastan_de_mas(concentradohogar2014),
                   gastan_de_mas(concentradohogar2016),
                   gastan_de_mas(concentradohogar2018))

over_sp <- ggplot()+
  geom_bar(aes(x=c("2012", "2014", "2016", "2018"), y=over_spenders), stat = "identity") +
  scale_y_continuous(labels = scales::percent) +
  xlab("Año") + ylab("")+
  labs(title="Porcentaje personas que gastan más de lo que ingresan")+
  theme_bw()

over_sp 

20% es un porcentaje para nada despreciable de gente que al parecer tiene gastos corrientes mayores que sus ingresos. Lo que es más, está métrica es muy consistente en los 4 años por lo que es justo asumir que estamos frente a un claro fenomeno y no simple error humano.

Datos del SAT

Aquí veremos una forma de aporvechar los datos que nos brinda el SAT para mejorar nuestro entendimiento de los ingresos de las personas. En específico nos interesa saber cómo se comportan los ingresos más extremos que reportan frente al SAT. Por supuesto, a la hora de modelar habrá que hacer las suposiciones adecuadas para poder compararlo con aquel ingreso corriente que se reporta en la ENIGH. De esto se habla en el resumen que acompaña a este código.

breaks <- c(1,2,3,4,5,6,8,10, 15,20,30,50,100,200,500,1000)*1000000


means <- c()
for (i in 1:(length(breaks)-1)) {
  means[i]<- mean(sat2012$I_DEC_TIAONCT1_AA[sat2012$I_DEC_TIAONCT1_AA >= breaks[i]
                                            & 
                                            sat2012$I_DEC_TIAONCT1_AA<breaks[i+1]], 
                  na.rm = T)
}

top_sat <- ggplot()+
  geom_line(aes(x=1:(length(breaks)-1), y=means), color = "light blue", size = 1) +
  scale_x_continuous(breaks = 1:(length(breaks)-1), labels = c("1-2 M", "2-3 M", "3-4 M", "4-5 M", "5-6 M", "6-8 M", "8-10 M", "10-15 M", "15-20 M", "20-30 M", "30-50 M", "50-100 M", "100-200 M", "200-500 M", "500-1000 M")) +
  scale_y_continuous(labels = scales::comma, breaks=seq(0,600000000,100000000)) +
  labs(title = "Ingreso medio segun ingresos reportados al SAT") +
  xlab("Ingreso medio") + ylab("Rango de ingreso")+
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

top_sat

Estimación Maxima verosimilitud Log-normal 2012

En seguida les presento una forma para estimar parámetros con y sin restricciones. La especificación de las restricciones es muy importante para lograr una solución al problema de optimización. En este caso, solo busco presentarles el flujo de trabajo y no parámetros 100% reales a los presentados en los artículos de Bustos.

mean_log <- function(x) {
  exp(x[1]+(x[2]^2)/2)
}

pll_log <- function(x){
  -sum(dlnorm(concentradohogar2012$ing_cor+.01, x[1], x[2], log = TRUE)*concentradohogar2012$factor)
}

par_log_NR<- optim(par = c(10, 1), pll_log, control=list(maxit=1000))$par


par_log_R <- solnp(pars = c(11, 2), # Valores inciales
               pll_log, # función a optimizar
               eqfun=mean_log, # función restricción
               eqB=92733.62,   # Valor a igualar
               LB=c(0,0), # Soporte mínimo para los parametros
               UB=c(20,5))$pars # Soporte máximo para los parametros
## 
## Iter: 1 fn: 373023071.8492    Pars:  10.53334  1.83813
## Iter: 2 fn: 369313786.9753    Pars:  10.47885  1.57183
## Iter: 3 fn: 367073387.1403    Pars:  10.55500  1.36961
## Iter: 4 fn: 366549763.1526    Pars:  10.60202  1.29591
## Iter: 5 fn: 366512324.9051    Pars:  10.61045  1.28615
## Iter: 6 fn: 366511843.7034    Pars:  10.61068  1.28593
## Iter: 7 fn: 366511843.4866    Pars:  10.61068  1.28593
## Iter: 8 fn: 366511843.4875    Pars:  10.61068  1.28593
## solnp--> Completed in 8 iterations
## Media muestral
weighted.mean(concentradohogar2012$ing_cor, concentradohogar2012$factor)
## [1] 37999.64
## Media No restringida
mean_log(par_log_NR)
## [1] 37254.33
## Media restringida
mean_log(par_log_R)
## [1] 92733.62

Estimación Maxima verosimilitud GB2 2012

Ahora observamos cómo se estima para otra distribución, en este caso la Beta generalizada de tipo II.

mean_gb2 <- function(x) {
  moment.gb2(1, x[1], x[2], x[3], x[4])
}

pll_gb2 <- function(x){
  -loglp.gb2(concentradohogar2012$ing_cor+.01, x[1], x[2],x[3], x[4], w=concentradohogar2012$factor)
}

par_gb2_NR<- optim(par = c(1.307071, 20594.790235, 2.341514, 1.843004), pll_gb2, method = "BFGS",control=list(maxit=1000))$par

par_gb2_R <- solnp(pars = c(1.5, 30000, 1, 1), # Valores inciales
               pll_gb2, # función a optimizar
               eqfun=mean_gb2, # función restricción
               eqB=92733.62,   # Valor a igualar
               LB=c(0,0,0,0), # Soporte mínimo para los parametros
               UB=c(5,100000,5,5))$pars # Soporte máximo para los parametros
## 
## Iter: 1 fn: 11.5445   Pars:      1.47087 30642.59739     1.32305     1.01227
## Iter: 2 fn: 11.5258   Pars:      1.48804 30640.46754     1.04509     0.94058
## Iter: 3 fn: 11.5258   Pars:      1.48801 30640.47383     1.04545     0.94067
## Iter: 4 fn: 11.5258   Pars:      1.48801 30640.47383     1.04545     0.94067
## solnp--> Completed in 4 iterations
## Media muestral
weighted.mean(concentradohogar2012$ing_cor, concentradohogar2012$factor)
## [1] 37999.64
## Media No restringida
mean_gb2(par_gb2_NR)
## [1] 38249.42
## Media restringida
mean_gb2(par_gb2_R)
## [1] 92733.62

Distribución empírica contra estimadas restringidas

Ahora nos gustaría saber qué tanta corrección hubo con respecto de la distribución original que nos arrojan los datos. Es decir, cómo afecta a las distribuciones estimadas el uso de una restricción en su media. De nuevo, una plática más extensa se encuentra en el resumen.

distribucion_empirica <- function(x, datos, factor){
  # x es el valor donde se quiere evaluar la función de distribución empírica
  # datos es el conjunto de datos en formato de vector que se quiere evaluar
  df <- data.frame(datos = datos, factor = factor, cumu = cumsum(factor))
  F_x <- c()
  for (i in 1:length(x)) {
    F_x[i] <- sum(df[df$datos<x[i], "factor"])/df[length(datos), "cumu"]
  }
  return(F_x)
}
distri_empi_restri <- ggplot() +
  geom_function(aes(colour = "Empirica"), 
                fun = distribucion_empirica, 
                args = list(datos = concentradohogar2012$ing_cor, 
                            factor = concentradohogar2012$factor), 
                size = 1) +
  geom_function(aes(colour = "GB2"), 
                fun = pgb2, 
                args = list(shape1 = par_gb2_R[1],
                            scale = par_gb2_R[2],
                            shape2 = par_gb2_R[3],
                            shape3 = par_gb2_R[4]),
                size = 1) +
  geom_function(aes(colour = "Log-Normal"), 
                fun = plnorm, 
                args = list(meanlog = par_log_R[1],
                            sdlog = par_log_R[2]), 
                size = 1) +
  scale_x_continuous(trans = 'log', limits = c(600, 4000000), breaks = 10^(seq(2,6,1)), labels = scales::comma) +
  scale_y_continuous(breaks = seq(0,1,.1)) +
  scale_colour_manual(name="Distribución",values=c("Empirica"="#3591d1", "GB2"="#f04546", "Log-Normal"="#f045d6")) +
  ylab("Probabilidad") + xlab("Ingreso (Escala log)") +
  labs(title = "Distribución empírica contra estimadas con restricción") +
  theme_bw() +
  theme(legend.position = c(.1,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))

distri_empi_restri

Restringida vs No restringida

Siguiendo el mimso proceso pero ahora presento las distribuciones con y sin restricciones. Vemos que las correciones son más o menos consistentes entre los modelos. Esto son buenas noticias pues podemos hacer algunas inferencias sin preocuparnos de casarnos con una distribución en específico.

p1 <- ggplot() +
  geom_function(aes(colour = "Log-Normal Restringida"), 
                fun = plnorm, 
                args = list(meanlog = par_log_R[1],
                            sdlog = par_log_R[2]), 
                size = 1) +
  geom_function(aes(colour = "Log-Normal No Restringida"), 
                fun = plnorm, 
                args = list(meanlog = par_log_NR[1],
                            sdlog = par_log_NR[2]), 
                size = 1) +
  scale_x_continuous(trans = 'log', limits = c(600, 4000000), breaks = 10^(seq(2,6,1)), labels = scales::comma) +
  scale_y_continuous(breaks = seq(0,1,.1)) +
  scale_colour_manual(name="Distribución",values=c("Log-Normal Restringida"="#f04546", "Log-Normal No Restringida"="#3591d1")) +
  ylab("Probabilidad") + xlab("Ingreso (Escala log)") +
  theme_bw() +
  theme(legend.position = c(.25,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))

p2 <- ggplot() +
  geom_function(aes(colour = "GB2 Restringida"), 
                fun = pgb2, 
                args = list(shape1 = par_gb2_R[1],
                            scale = par_gb2_R[2],
                            shape2 = par_gb2_R[3],
                            shape3 = par_gb2_R[4]),
                size = 1) +
  geom_function(aes(colour = "GB2 No Restringida"), 
                fun = pgb2, 
                args = list(shape1 = par_gb2_NR[1],
                            scale = par_gb2_NR[2],
                            shape2 = par_gb2_NR[3],
                            shape3 = par_gb2_NR[4]),
                size = 1) +
  scale_x_continuous(trans = 'log', limits = c(600, 4000000), breaks = 10^(seq(2,6,1)), labels = scales::comma) +
  scale_y_continuous(breaks = seq(0,1,.1)) +
  scale_colour_manual(name="Distribución",values=c("GB2 Restringida"="#f04546", "GB2 No Restringida"="#3591d1")) +
  ylab("Probabilidad") + xlab("Ingreso (Escala log)") +
  theme_bw() +
  theme(legend.position = c(.2,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))

grid.arrange(p1,p2,ncol=2)

Ingreso acumulado por Decil empírico contra estimados

Por último, revisemos una muestra rápida de cómo se ve la desigualdad estimada cuando asumimos que los modelos antes descritos son los correctos frente a la desigualdad que existe entre la muestra del ENIGH.

## Calculamos el acumulado por decil mediante simulación
x <- rlnorm(100000, par_log_R[1], par_log_R[2])
deciles <- quantile(x, probs = seq(0,1,.1))
total<-sum(x)

medd <- c()
for (i in 1:10) {
  medd[i] <- sum(x[x > deciles[i] & x <= deciles[i+1]])/total
}

## Construcción de la tabla para alimentar el ggplot
n_percentil <- 10
df <- tibble("ENIGH 2012" = ingreso_acumulado_percentil(concentradohogar2012, n_percentil),
             "Log-Normal" = medd,
             index = fct_rev(as.factor(1:n_percentil))) 

df <- df %>% gather(key = "Fuente", value = "Ingreso", -index)

## Construcciónd de la gráfica
ingreso_acu_decil_ENIGH_Log <- ggplot(df, aes(fill = index, x = Fuente, y = Ingreso)) +
  geom_bar(position="stack", stat="identity", width=0.7) +
  ylab("Porcentaje de Ingreso acumulado") + xlab("Modelo") +
  scale_fill_viridis(discrete = T, name = "Decil") +
  scale_y_continuous(labels = scales::percent) +
  theme_bw() +
  theme(legend.background = element_rect(fill = "transparent", colour = "transparent"))

ingreso_acu_decil_ENIGH_Log

Podemos observar que la desigualdad se acrecenta considerablemente. Esto por que el decil más rico acapara un porcentaje aún mayor de el ingreso corriente.

Simulación nube de gas

Se generan los individuos y se especifica el total del dinero en la simulación

N = 500
din_total = 5*10^5

Se hace un dataframe para guardar la información de las simulaciones

dataind = data.frame(
  "index" = c(1:N),
  "dinero" = din_total/N
)

Se realiza la simulación

for (i in 1:(4*10000)){
  # elegir los dos agentes de la transacci?n
  primer_agente = round(runif(1, 1, 500))
  segundo_agente = round(runif(1, 1, 500))
  while (primer_agente == segundo_agente){
    segundo_agente = round(runif(1, 1, 500))
  }
  # determinar el ganador
  ganador = rbinom(1, 1, 0.5)
  
  if (ganador == 1){
    indice_ganador = segundo_agente
    indice_perdedor = primer_agente
  } else {
    indice_ganador = primer_agente
    indice_perdedor = segundo_agente
  }
  
  # calcular el dinero a ser transferido usando la regla tres
  transaccion = runif(1, 0, 1) * din_total/N
  
  # checar si el perdedor tiene suficiente dinero para que se
  # realice la transacci?n
  
  if (dataind[indice_perdedor, "dinero"] < transaccion){
    next
  }
  
  # llevar a cabo la transacci?n
  
  dataind[indice_ganador, "dinero"] = dataind[indice_ganador, "dinero"] + transaccion
  dataind[indice_perdedor, "dinero"] = dataind[indice_perdedor, "dinero"] - transaccion
}

Una vez que se ha realizado la simulación, se grafica un histograma para obtener la densidad de la distribución del ingreso para los agentes y sobre este histograma, se grafica la función de densidad propuesta para modelar la distribución de este ingreso

hist(dataind[,"dinero"],freq = FALSE)
lines(seq(1:7000),  (densidad = (1/1000) * exp(-seq(1:7000) / 1000)))